      SUBROUTINE INTER(FI,DFI,CX,CY,X,Y,CINC,SOL,DSOL,KODE,INC,NNINC,
     &                 NLD,NG,NPI,NN,N0,NINCL,C0)

C ESTA SUBROTINA CALCULA O VALOR DO POTENCIAL PARA PONTOS INTERNOS


      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'param_dim.inc'

!       COMMON N_Eq,NN,L,LEC,IMP,nu_msc,nu_msd,nu_trip,nu_trif

      !       COMMON /incdata/ C0,N0,NINCL
      COMMON /PG/ GI(8,2),OME(8,2)

      DIMENSION FI(NX),DFI(NX),CX(NPIX),CY(NPIX),X(NX),Y(NX)
      DIMENSION SOL(NPIX),DSOL(NPIX,2)
      DIMENSION KODE(NX),INC(NX,2)

      DIMENSION NNINC(NINC_max),CINC(NINC_max),NLD(NINC_max+1,2)

! ! ! C
! ! ! C REORDENAR OS VETORES FI E DFI: 
! ! ! C                            POTENCIAL EM FI 
! ! ! C                            DERIVADAS EM DFI
! ! !       DO 20 I=1,N0
! ! !          IF(KODE(I).GT.0) THEN 
! ! !            CH=FI(I)
! ! !            FI(I)=DFI(I)
! ! !            DFI(I)=CH
! ! !          END IF
! ! !    20 CONTINUE
! ! ! 
! ! ! C     Copia os valores calculados do 
! ! ! C     potencial nas interfaces para FI
! ! ! 
! ! !          nna = nn
! ! !          nna0= n0
! ! ! 
! ! !       do iinc =1,nincl
! ! !          nni = nninc(iinc)
! ! !          
! ! !          do i=1,nni
! ! !             fi(nna0+i)= dfi(nna0+i)
! ! !             dfi(nna0+i)= dfi(nna+i)
! ! !         enddo
! ! ! 
! ! !          nna = nna  + nni
! ! !          nna0= nna0 + nni
! ! !       enddo
C
C CALCULO DO VALOR DO POTENCIAL PARA PONTOS INTERNOS
C

      ALFA = 100.0
      B = 0.0
      BJ = 0.0
      db1 = 0.
      db2 = 0.


C     *************************************     
C     *************************************     

C     Implementacao  para inclusoes em desenvolvimento ...

C     *************************************     
C     *************************************     


      if(npi.gt.0) then

C     Flag para indicar o calculo de resultados em pontos internos.
         IPTI=1

         do 50 idom= 1,nincl+1

            if (idom.gt.1) write(*,*) 'idom,c,sig',
     &                       idom,cinc(idom-1),-1.0*c0/cinc(idom-1)

C     Armazena a condutividade do meio em cond
            if(idom.eq.1) then
               cond= c0
            else
               cond= cinc(idom-1)
            endif

C     Determina os nos inicial e final dos contornos de cada dom�nio
            IF(IDOM.EQ.1) THEN
               NPTI_I= nld(idom,1) - nn
               NPTI_F= nld(idom,2)- nn
               
               NEL_I= 1
               NEL_F= NN
            ELSE
               NPTI_I= nld(idom,1) - nn
               NPTI_F= nld(idom,2) - nn

               IF(IDOM.EQ.2) NEL_F = N0
               NEL_I= NEL_F + 1
               NEL_F= NEL_F + NNINC(IDOM-1)

            ENDIF

!             write(*,*) 'Elemento Inicial e Final do Dom.',idom
!             write(*,*)nel_i,nel_f,NPTI_I,NPTI_f

		write(*,*)'npti_i,npti_f:',npti_i,npti_f
            DO 40 K=npti_i,npti_f
C        Zera os vetores que irao acumular as repostas nos pontos internos.



               SOL(K)=0.0
               DSOL(K,1)= 0.0
               DSOL(K,2)= 0.0
C     Faz um Loop sobre os elementos.
               DO 30 J=nel_i,NEl_f

                  NI= INC(J,1)
                  NF= INC(J,2)
                  sig= 1.d0
C     Inverte as incidencias para PTOS. INT. nas inclusoes
                  if(idom.gt.1)then
                     ICH= NI
                     NI=  NF
                     NF= ICH
                     sig= -1.0*c0/cinc(idom-1)
                  endif
                  CALL INTE(CX(K),CY(K),X(NI),Y(NI),X(NF),Y(NF),
     &             HJ,DHJ1,DHJ2,GJ,DGJ1,DGJ2,B,BJ,db1,db2,alfa,IPTI,NG)     
!                   CALL INTE(CX(K),CY(K),X(NI),Y(NI),X(NF),Y(NF),
!      &                      HJ,DHJ1,DHJ2,GJ,DGJ1,DGJ2,IPTI,NG)
                  SOL(K)=    SOL(K) + sig*DFI(J)*GJ - FI(J)*HJ 
                  DSOL(K,1)=    DSOL(K,1) + sig*DFI(J)*DGJ1 - FI(J)*DHJ1 
                  DSOL(K,2)=    DSOL(K,2) + sig*DFI(J)*DGJ2 - FI(J)*DHJ2 
            
 30            CONTINUE

               
               SOL(K)=    0.5* SOL(K)  /(3.1415926)

C     Calcula as componentes do fluxo como - cond * Grad(u)
               DSOL(K,1)= -1.*cond*0.5*DSOL(K,1)/(3.1415926)
               DSOL(K,2)= -1.*cond*0.5*DSOL(K,2)/(3.1415926)
                        
 40         CONTINUE
 50      continue
      endif
      
      RETURN
      END
